The statistics of diffusive flux 
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We calculate the explicit probability distribution function for the flux between 
sites a in a simple discrete-time diffusive system composed by independent random 
walkers. We highlight some of the features of the distribution and we discuss its 
relation to the local instantaneous entropy production in the system. Our results 
are applicable both to equilibrium and non equilibrium steady states as well as for 
certain time dependent situations. 

PACS numbers: 



Random walkers are used extensively as phenomenological models for diffusive 
processes Q|. The connection between physical diffusion and random walks arises from the 
fact that the concentration in a system containing many identical independent random walk- 
ers satisfies Pick's law. That is, the mean flux between sites is proportional to the difference 
in the mean number of walkers at each site, plus, perhaps, a convective term if the walks 
are not symmetric. However, to the best of our knowledge, the statistics of the flux in such 
a system have not been characterized (beyond the mean, of course). This is not the case 
in more complicated systems, for which a large amount of work on the statistics of currents 
has been carried out. Recent examples include statistics and large deviation theory for the 
fluctuations of the current in lattice gases ; the integrated current distribution for 

the steady-state of the one-dimensional zero-range process p]; the joint probability function 
for the occupation number and the current through the system in the asymmetric simple 
exclusion process with open boundaries 0] , to mention but a few. 

In addition to being an alternative approach to the description of transport properties 
of simple diffusive systems, the statistics of flux are of relevance in many diffusion limited 
processes. For example, in a diffusion reaction system in which the reactants are initially 
separated, fluctuations in the flux of reactants into the interfacial region give rise to fiuctu- 



ations in the position of the reaction front and in the net reaction kinetics 

Furthermore, viewed from the thermodynamic side, there has been a great deal of work 
relatin g th e statistics of flux to the entropy production of non equilibrium steady states 
|llLll2lll3l|. From this perspective, diffusive systems modeled by random walkers are among 
the simplest examples that can be used to test and extend such ideas. 

In this work we obtain an explicit expression for the probability distribution of the flux, 
that is, the net number of particles that hop between neighboring sites in a given time 
step, in a system containing independent discrete time random walkers. This system has 
the added bonus that the results are applicable to both equilibrium and non-equilibrium 
situations, where the latter can be achieved by imposing concentration differences on the 
boundaries or by subjecting the system to an external field (by considering biased random 
walkers) or both. 

For definiteness, in this work we consider a discrete one dimensional system in which 
independent discrete time random walkers evolve synchronously according to the following 
rules: at each time step, a walker moves to the site on its right or to the site on its left with 
probabilities p and q respectively, or stays at its site with probability r = 1—p — q. We note 
that the discreteness in time allows us to consider fluxes that arise from the simultaneous 
hopping of many particles. Such multiparticle events are not present in continuous time 
versions of the system; we discuss the continuous limit of the problem further on. 

The system in which the random walkers evolve consists of / + 2 sites and we impose as 
boundary conditions that the number of walkers at sites i = and i = / + 1 be Poisson 
distributed with fixed mean values Cq and q+i respectively [l^ . 

To compute the probability distribution of particle flux between neighboring sites, labeled 
i and z + 1, we require J+ and J_, the number of particles that jump from site i to site i + 1 
and the number of particles that jump from i + 1 to i respectively. In terms of these, the 
total flux J between these sites will be given by J = J+ — J_. 

We denote by m and n the ocupancy of sites i and i + 1 respectively. Then, the probability 
that the net flux between these sites is J can be expressed as 



P 



(J) = ^p{J\m,n)pi^i+i{m,n), (1) 



where pj^j+i(m, n) is the probability of finding exactly m and n particles at sites i and i + 
and p{J\m,n) is the probability of having a total flux J between these sites given those 
occupancy numbers. 
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Since the walkers are independent, we have 



oo 



p(J\m,n) = p+(J + J_|m)p_( J_|ra), (2) 

J-=0 



where 

(TTl \ 
J l/Mi-pr~'+ (3) 

and 

P-{J-\n)= i^^q'-{l-qr-'- (4) 

are the probabihties that J+ particles jump to the right from a site containing m particles 
and J_ particles jump to the left from a site containing n particles. 

In what follows, it will prove useful to work with the generating function of expression 



(j21), defined as p(z|m, n) := X]jL-oo z^'p{.J\''^^ which is given by 



p{z\m,n) = (l - q+^^ (l-p + zp)"". (5) 

Next, we note that this system has the remarkable property that the joint probability 
distribution_»j j+i(m, ra) of the occupancy numbers factorizes into a product of Poisson dis- 



tributions 



15| . This property is analogous to the factorization of the steady state occupancy 



distribution which occurs in certain zero range processes under similar boundary conditions 
jlfi| . Thus, the joint occupancy distribution pj^j+i(m, can be written as 

Pi^i+i{m,n) = Lpi{m)(fi+i{n), (6) 

where 

= ^e--cf , (7) 

and Cj is the mean number of particles at site i. Actually, under certain circumstances, 
namely initial conditions already characterized by independent Poisson distributions, the 
joint occupancy distribution can be expressed as the product of Poisson distributions 
throughout the evolution of the system. In this case, the mean occupation numbers can 
be obtained as the solution to the discrete diffusion equation 



Ci{t + 1) - Ci{t) = pci_i{t) + qci+i{t) - (p + q)ci{t) 



(8) 
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with appropriate initial and boundary conditions. However, just as for zero range processes 
jlG^ . the system always reaches a unique factorizable steady state, independently of the 
initial conditions, for which the parameters q are given by 

i+i 



1 



i+i 

V -1 



Co 



- Q+i + (q+1 - Co) 



(9) 



which is the steady state solution of equation (jH)). 

Thus, we are in a position to evaluate the generating function Pi{z) of the distribution of 
fluxes between sites i and i + 1: 



(10) 



From this expression, the moments and the cumulants of the distribution may be readily 
evaluated. For example, all the odd cumulants are given by 



i^2n-i = pci - qci+i, n = 1, 2, 



(11) 



whereas the even cumulants are 



i^2n = pCi + qci+i, n = l,2,. 



(12) 



In particular, using the explicit expression for the concentration profile for the steady states 
given in equation 0, we find 

1 



pi+i _ qi 



-[(p-g)(V+^-QW+')] 



and 



[{P + g)(cop'+' - Q+ig'+i) + 2m(q+i - c^)p'q'-'] 



(13) 



(14) 



Furthermore, using the explicit expressions for the moments we can calculate the skewness 
7i and excess (or kurtosis) 72 of the distribution, which are given by 

{{J-{J)f) (pc,-gc,+i) 



71 



0-^ 



{pCi + gci+i)2 



and 



72 = 3 



a'* 



(15) 



(16) 
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Thus, only in the hmit pcj + gcj+i ^ 1 is the distribution of flux essentially Gaussian. In fact, 
the expression for Pi{z) can be inverted to obtain an explicit expression for Pi{J). Rewriting 



-{2,/pqCiCi+i 




(17) 



and comparing with the generating function of the modifled Bessel functions 

exp ia; (t + ij = ^ Is{x)f, 



yields 



-pci-qci+i 



PCj 
qCi+1 



Ij{2y/pqCiC 



i+l) 



(18) 



(19) 



This distribution, shown in Figd shares some of the properties of the Gaussian distribu- 
tion: It is characterized by only two parameters, say pci and gQ+i in this case. It satisfles 
a generalized stability condition, in the sense that the distribution of the sum of indepen- 
dent (integer) random variables distributed according to Pi{j), has the same form as Pi{j) 
with appropriately rescaled parameters. That this is the case is obvious from the explicit 
expression of the generating function, but it could have been foreseen from the actual pro- 
cess we are describing. Indeed, the complete derivation presented above applies also for the 
distribution of flux between adjacent cell in systems of arbitrary dimension. For these, the 
one dimensional system we are considering corresponds to a projection onto one of the axes. 
Thus, the flux distribution we obtain can itself be thought of as that of the sum of many 
independent fluxes with distributions similar to Pi{j)- 

The large \J\ behavoir of Pi{J), which is where it most clearly differs from the Gaussian, 
is easily amenable to evaluation. From the series expansion of the Bessel functions we 
immediately flnd that 



0_ _j_ PgCjCj+l _|_ . . . ) 



p,{J) ~ < 



^ \J\\ \J\ ^ 



when J — > oo 



when J — > — oo. 



(20) 



Another limiting behavior of Pi{J) worth mentioning is that corresponding to the con- 
tinuous time limit. This limit is achieved by assuming that the hopping probabilities can 
be expressed as p = a5t and q = (36t, where 6t is the time interval between succesive steps. 
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FIG. 1: Semilogarithmic plots of Pi{J) with variance equal to 2 and a) (J) = 1; b) (J) = (dot- 
dashed lines). Gaussian distributions with the same variance and mean as Pi{j) in each case are 
plotted for comparison. 



Then, keeping only the terms up to linear order in 6t we find 



aCiSt J = 1 

1 - {aci + (3ci+i)6t J = 
(3ci+i5t J = —1. 



(21) 



All other values of J have probabilities of higher order in 6t and are, therefore, negligible 
as 5t — > 0. The above expression merely reflects the fact that in the continuum limit, the 
transport of particles amongst neighboring sites during a time interval 6t is, at most, a single 
particle process. The statistical parameters in this limit can be calculated directly from the 
expressions obtained above; thus, for example, the mean number of particles that flow from 
site i to site i + 1 in a time interval 6t is (ac, — /?Ci+i)5t as was to be expected. 

A further characteristic that Pi{J) shares with the Gaussian distribution is that the ratio 
Pi{J)/pi{—J) is a pure exponential. This quantity is of interest because it can be associated 
to entropy production of non-equilibrium steady states. Indeed, given the explicit form of 
Pi{J) between sites i and i -|- 1 in a system, and the symmetry properties of the Bessel 
functions, we have: 



In 



J In 



PCj 

qci+i 



(22) 
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The average of this quantity is 

— (23) 

which has the famihar non-negative form that appears in the i!f-theorem, and that we argue 

can be considered as an instantaneous local entropy production for this system. It should 
)e remarked that this is not the same quantity that appears in the fluctuation theorems 
111 Q| , which are concerned with the large fluctuations of the time- integrated flux. 
Also, from the definition of Si in equation (j22j) it is again apparent that this quantity is a 

random variable, the statistics of which are simply related to the statistics of flux described 

above. Thus, for example, the variance of Si will be given by 

2 



(Si) - {Si^ = (pci + qci+i) 



In 



pCi 



(24) 



To justify our identification of (sj) as an average local entropy production, first note that 
the equilibrium for this system is defined as the steady state with zero average flux. In this 
situation, the concentration profile will be given by 

c^ = Co{p/qy, (25) 

where Cq is the concentration imposed on the boundary site i = 0. On the other hand, 
from the thermodynamic perspective the system corresponds to noninteracting particles in 
an external field h. For such a system the concentration profile is given by the barometric 
equation 

c. = ef^^^-'^l (26) 

where /z is the chemical potential and /3 is the inverse temperature. Comparing both ex- 
pressions for the concentration profiles at equilibrium leads to the identification 

/?// = lnco; (3h = \n-. (27) 

q 

Of course, in general the system is not in equilibrium; however, since the occupancy is 
characterized by independent Poisson distributions, it can be considered as being in local 
equilibrium. Thus, we can rewrite equation ()23|) as: 



(Sr) = {J)^ 



P 

In Cj — In Cj+i + In - 
1. 



-(J),[A,/5(/x, + $(0)], (28) 



where $(«) = —ih is a "potential energy" and Aj is the difference operator. Equation 
is a discrete analog of the local entropy production of linear thermodynamics. 
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At this point it is worth emphasizing that the average entropy production proposed in 
equation (j^^ is obtained from a single step local statistic. Thus, this quantity should 
be interpreted as the entropy produced at that time step, averaged over an ensemble of 
statistically identical systems, which is the usual interpretation of statistical mechanics. 
This is especially relevant in the time dependent case, for which all of the above results also 
hold, as mentioned above, given that the initially the sites of each system of the ensemble 
are occupied according to independent Poisson distributions. For these it makes no sense 
to study time averaged quantities of single systems and relate the results to the ensemble 
averages. On the contrary, in the steady state the ensemble average of Sj, Eq. (P^ . will 
coincide with the time average of the flux in a single system, yielding the connection usually 
assumed by the ergodic hypothesis. This will not be true for higher moments of Sj due to 
correlations in the flux at different times. 

If we restrict ourselves to the steady states, the average flux is constant and we can 
calculate the total entropy production in the system as 

I I / c \ 
Stotai = y~'(-5)i = y](pci - qci+i) In ( — I (29) 

= pi+i ! qi+i i^P - l^^'^^P'^' - ^'+1^'^')] [(^ + 1) In + 111 (co/q+i)] , 

where Cq and q+i are the concentrations imposed on the boundaries of the system and we 
used the explicit expression for the steady state flux given in equation 0. Inparticular, if 
we take the limit q = p, the above expression coincides with that obtained in [15[ , where the 
entropy production is interpreted in terms of a time asymmetry in the dyna mical randomness 
between the forward and backward paths of the diffusion process 

HQ. 

Finally, the limit in which q —>■ say, is worth discussing. This totally asymmetric hopping 
case is again of some importance in the context of zero range processes. The distribution 
of flux in this case is easy to evaluate from the explicit expression in equation (fT^ . from 
which we obtain Pi{J) = for J < and Pi{J) is Poissonian with mean pci for J > 0. 
Further, as discussed above, a field can be associated to the ln(p/g), which gives rise to a 
"potential energy" ; thus, in the limit when q vanishes, motion against the field is impossible 
and motion along the field "dissipates" an infinite amount of energy. In accordance, the 
mean entropy production diverges in this limit. 

In summary, we have obtained the explicit probability distribution function of single step 
flux between adjacent sites in a diffusive system composed of independent discrete random 
walkers. The key to solving the problem is the fact that the joint site occupancy distribution 
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for this system factorizes into independent Poisson distributions even in non equilibrium 
situations. We also obtain the statistical parameters that characterize the distribution and 
discuss some of its limiting behaviors. In addition, we use the distribution to evaluate a 
single step local entropy production, whose average can be related to the usual expressions 
from nonequilibrium thermodynamics. In this context, an interesting extension of this work 
would be to introduce an additional conserved quantity which can be distributed among the 
random walkers at each site. This energy-like quantity may or may not affect the jumping 
rates (although from a physical perspective it probably should), and would allow the study 
of coupled transport phenomena in these extremely simple settings. We have thus far been 
unable to pose such scenario in a tractable manner. A different avenue of research would 
be to determine multiple time and site flux statistics in systems containing many random 
walkers, characterizing, for example, the correlations in flux at two sites of the system as a 
function of time. 

This work was partially supported by grant IN-100803 of DGAPA UNAM. We thank D. 
Sanders, F. Leyvraz and M. Aldana for the many useful comments and suggestions on this 
manuscript. 

APPENDIX 

The key result that permits the explicit calculation of Pi{J) is that the joint occupation 
probability distribution of the system factorizes into Poisson distributions. For the sake of 
completeness, we show how this comes about. 

We require the determination of the joint occupation probability distribution, that is, the 
probability pt{no, rii, . . . , n^+i) of flnding no, ni, . . . , n;+i particles, at the sites 0, 1, ...,/ + 1 
respectively at time t, with the boundary conditions described in the text. The evolution 
equation for this quantity is 

Pt+i{no,ni,...,ni+i) = e^'^'^-^e-'^'+i-^ ( pt{mo,mi, . . . ,mi+i) (30) 

no! n;+i! J 

X l[6{m + lU + ll^i - {It + ni) 
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where 

oo oooooo oooooooo oo 

/ = E' E E E -E E EE ' E 

mo=0 mi=0 mi+i=0 I- =0 I- =0 1-^^=0 1+ =0 1+ =0 11+1=0 

and we recall that r = 1 — p — q. Equation ()3Hl looks messy, but it is actually very easy to 
understand: if we consider the LHS of this equation without the delta functions, we note 
that it contains all the possible events from all the possible configurations at time t, weighted 
by the probability that the events occur. The deltas restrict this immense sum to the terms 
which end up in the occupancy configuration specified on the RHS of the equation. 

In expression (jHlj) we have already imposed Poisson occupancy distributions with con- 
centrations Co and Cj+i on the sites and / + 1 as boundary conditions for the system. This 
will enable us to drive the system into "nonequilibrium" steady states by imposing different 
concentrations on the boundary, or considering biased random walks (i.e. p ^ q), or both. 

This evolution equation can be dealt with by evaluating the multivariate generating 
function of the joint distribution, namely: 

oo oo oo 

Pt{zo, . . . , zi+i) = X] X] • • • ^o°^i' • • • ^z+i'-P*(^o, ni, . . . , rii+i). (31) 

no=Oni=0 n;+i=0 

After some simple but tedious algebra, which consists in regrouping terms and adding up 
binomial expansions, the generating function can be shown to satisfy the recursion relation 

Pt+i(zo, . . . , zi+i) = e-^«(i-^«)e-'='+^(i-^'+^) (32) 

oo oo 

X J^... ^ pt{mo,mi,...,mi+i) 

mo=0 m,+i=0 

X [1 + (z, - l)p]™" [q + Z2P + ZirY"' 

X [zi.iq + P + zirr^ [1 + {zi - l)g]"^'+^ 
i-i 

X \\[zi+ip + Zi-iq + ZirY'\ 

1=2 

Comparing with the definition of the generating function, we note that this expression is 
equivalent to 

p,+i(^o, ^1, . . . , ^z, zi^i) = e-'=o(i-^o)e-'='+^(i-^'+^) X (33) 
Pt (l + {zi - l)p, Z2P + q + zir, zsp + Ziq + • • • 

. . . , zip + zi_2q + zi_ir, zi-iq +p + zir, 1 + {zi - l)q). 
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To solve this equation we make the ansatz that the distribution is given by the product 
of Poisson distributions, and thus, its generating function can be written as 

Mzo,...,zi+,) = e-^'^-oc^^'^^'-^'\ (34) 

Substitution of this expression into Eq. ljH^ leads to the following equations for the mean 
occupation numbers Ci{t) which characterize the Poisson distributions at each site: 
l+l 

Ci{t + 1)[1 - Zi] = Co[l - ^o] + Q+i[l - zi+i] + Co(t)p[l - zi] + ci(t)[l - {Z2P + q + Zir)] 

i=0 

l-l 

+ ^Ci{t)[l - {Zi+ip + Zi_iq + Zir)] + ci{t)[l - {zi_iq + p + Zir)] + Q+i(t)g[l - Zi] 

i=2 

If we now impose the boundary conditions co(t) = cq and Q+i(t) = q+i for all times t, the 
above equation may be rewritten as 
I I 

^Ci{t + 1) - J^Ci(t) = {cop - ci(t)g) - {ci{t)p - Q+ig) (35) 

i=l i=l 
I 

+ X^[ci(i + 1) - Ci-i(t)p - Ci+i{t)q - Ci{t) + Ci{t){p + q)]zi 
1=1 

which must be valid for all values of the variables z^. This is achieved if, in addition to the 
boundary conditions, the parameters Ci{t) satisfy the discrete diffusion equation 

Ci(t + 1) - Ci(t) = pci-i(t) + qci+i{t) -{p + q)ci(t). (36) 

Thus, if we have an ensemble of systems in which each site is initially occupied according 
to independent Poisson distributions, then the occupancy distribution will remain being the 
product of Poisson distributions throughout the evolution of the process. 

Furthermore, independently of the initial conditions, the system will reach a unique 
stationary state characterized by the product distribution with the parameters q solutions 
of the stationary state of equation (jH^ : 



1 



l+l 



Co ( - Q+i + (q+1 - Co) 



(37) 



- 1 

With this expression, we can calculate the explicit values of the mean flux and its variance 
as functions of p, q, the imposed boundary concentrations Cq and q+i, and, for the variance, 
the position along the system: 

(Jn) = ^p^T^i [iP - - ^'+1^'^')] (38) 
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2 (P + Q) , J \ , ^ f Ci+i - Co \ i , . 

Finally, it is worth mentioning that the above derivation can also be extended to include 
site dependent jump rates. 
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